#include <misc.h>
#include <preproc.h>

module surfrdMod

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: surfrdMod
!
! !DESCRIPTION:
! Contains methods for reading in surface data file and determining
! two-dimensional subgrid weights as well as writing out new surface
! dataset. When reading in the surface dataset, determines array
! which sets the PFT for each of the [maxpatch] patches and
! array which sets the relative abundance of the PFT.
! Also fills in the PFTs for vegetated portion of each grid cell.
! Fractional areas for these points pertain to "vegetated"
! area not to total grid area. Need to adjust them for fraction of grid
! that is vegetated. Also fills in urban, lake, wetland, and glacier patches.
!
! !USES:
  use shr_kind_mod, only : r8 => shr_kind_r8
  use abortutils  , only : endrun
  use clm_varpar  , only : lsmlon, lsmlat
  use clm_varpar  , only : nlevsoi, numpft, &
                           maxpatch_pft, maxpatch_cft, maxpatch, &
                           npatch_urban, npatch_lake, npatch_wet, npatch_glacier
  use clm_varsur  , only : wtxy, vegxy
  use clm_varsur  , only : pctspec
  use ncdio
  use clmtype
  use spmdMod                         
  use clm_varctl,   only : scmlat, scmlon, single_column
  use decompMod   , only : get_proc_bounds,gsMap_lnd_gdc2glo,perm_lnd_gdc2glo
!
! !PUBLIC TYPES:
  implicit none
  save
!
! !PUBLIC MEMBER FUNCTIONS:
  public :: surfrd  ! Read surface dataset and determine subgrid weights
  public :: surfrd_get_grid  ! Read surface dataset into domain
  public :: surfrd_get_latlon  ! Read surface dataset into domain
  public :: surfrd_get_frac  ! Read land fraction into domain
  public :: surfrd_get_topo  ! Read topography into domain

!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
! Updated by T Craig
!
!EOP
!
! !PRIVATE MEMBER FUNCTIONS:
  private :: surfrd_wtxy_special
  private :: surfrd_wtxy_veg_rank
  private :: surfrd_wtxy_veg_all
  private :: surfrd_wtxy_veg_dgvm
  private :: surfrd_mkrank
!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd
!
! !INTERFACE:
!  subroutine surfrd(vegxy, wtxy, lfsurdat, domain)
  subroutine surfrd(lfsurdat, domain)
!
! !DESCRIPTION:
! Read the surface dataset and create subgrid weights.
! The model's surface dataset recognizes 5 basic land cover types within
! a grid cell: lake, wetland, urban, glacier, and vegetated. The vegetated
! portion of the grid cell is comprised of up to [maxpatch_pft] PFTs. These
! subgrid patches are read in explicitly for each grid cell. This is in
! contrast to LSMv1, where the PFTs were built implicitly from biome types.
!    o real edges of grid
!    o integer  number of longitudes per latitude
!    o real latitude  of grid cell (degrees)
!    o real longitude of grid cell (degrees)
!    o integer surface type: 0 = ocean or 1 = land
!    o integer soil color (1 to 20) for use with soil albedos
!    o real soil texture, %sand, for thermal and hydraulic properties
!    o real soil texture, %clay, for thermal and hydraulic properties
!    o real % of cell covered by lake    for use as subgrid patch
!    o real % of cell covered by wetland for use as subgrid patch
!    o real % of cell that is urban      for use as subgrid patch
!    o real % of cell that is glacier    for use as subgrid patch
!    o integer PFTs
!    o real % abundance PFTs (as a percent of vegetated area)
!
! !USES:
    use clm_varctl  , only : allocate_all_vegpfts
    use pftvarcon   , only : noveg
    use fileutils   , only : getfil
    use domainMod , only : domain_type
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
!    integer , intent(out) :: vegxy(:,:)   ! PFT
!    real(r8), intent(out) :: wtxy(:,:)  ! subgrid weights
    character(len=*), intent(in) :: lfsurdat               ! surf filename
    type(domain_type),intent(in) :: domain ! domain associated with wtxy
!
! !CALLED FROM:
! subroutine initialize in module initializeMod
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein, Sam Levis and Gordon Bonan
!
!EOP
!
! !LOCAL VARIABLES:
    character(len=256) :: locfn                          ! local file name
    integer  :: ncid,dimid,varid                         ! netCDF id's
    integer  :: begg,endg   
    logical  :: found                                    ! temporary for error check
    integer  :: iindx, jindx                             ! temporary for error check
    integer  :: ier                                      ! error status
    character(len=32) :: subname = 'surfrd'              ! subroutine name
!-----------------------------------------------------------------------

    call get_proc_bounds(begg,endg)
    allocate(pctspec(begg:endg))

    vegxy(:,:)   = noveg
    wtxy(:,:)  = 0._r8
    pctspec(:) = 0._r8

    if (masterproc) then
       write (6,*) 'Attempting to read surface boundary data .....'
       if (lfsurdat == ' ') then
          write(6,*)'lfsurdat must be specified'; call endrun()
       endif
    endif

    ! Read surface data

    if (masterproc) then
       call getfil( lfsurdat, locfn, 0 )
       call check_ret( nf_open(locfn, 0, ncid), subname )
    end if

    ! Obtain surface dataset special landunit info

    call surfrd_wtxy_special(ncid, domain)

    ! Obtain surface dataset vegetated landunit info

#if (! defined DGVM)     
    if (allocate_all_vegpfts) then
       call surfrd_wtxy_veg_all(ncid, domain)
    else
       call surfrd_wtxy_veg_rank(ncid, domain)
    end if
#else
    call surfrd_wtxy_veg_dgvm(domain)
#endif

    if ( masterproc )then
       call check_ret(nf_close(ncid), subname)
       write (6,*) 'Successfully read surface boundary data'
       write (6,*)
    end if

  end subroutine surfrd

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_get_grid
!
! !INTERFACE:
  subroutine surfrd_get_grid(domain,filename)
!
! !DESCRIPTION:
! Read the surface dataset grid related information:
! o real edges of grid
! o integer  number of longitudes per latitude
! o real latitude  of grid cell (degrees)
! o real longitude of grid cell (degrees)
! For offline mode only the the grid does not have to be global.
! If grid is read in from dataset, grid is assumed to be global but
! does not have to be regular.
! If grid is generated by model, grid does not have to be global but 
! must then define the north, east, south, and west edges:
!
! o edges(1)    = northern edge of grid (degrees): >  -90 and <= 90
! o edges(2)    = eastern edge of grid (degrees) : see following notes
! o edges(3)    = southern edge of grid (degrees): >= -90 and <  90
! o edges(4)    = western edge of grid (degrees) : see following notes
!
!   For partial grids, northern and southern edges are any latitude
!   between 90 (North Pole) and -90 (South Pole). Western and eastern
!   edges are any longitude between -180 and 180, with longitudes
!   west of Greenwich negative. That is, western edge >= -180 and < 180;
!   eastern edge > western edge and <= 180.
!
!   For global grids, northern and southern edges are 90 (North Pole)
!   and -90 (South Pole). The western and eastern edges depend on
!   whether the grid starts at Dateline or Greenwich. Regardless,
!   these edges must span 360 degrees. Examples:
!
!                              West edge    East edge
!                            --------------------------------------------------
!  (1) Dateline            :        -180 to 180       (negative W of Greenwich)
!  (2) Greenwich (centered):    0 - dx/2 to 360 - dx/2
!
!    Grid 1 is the grid for offline mode
!    Grid 2 is the grid for cam and csm mode since the NCAR CAM
!    starts at Greenwich, centered on Greenwich
!
! !USES:
    use clm_varcon, only : spval
    use domainMod , only : domain_type,domain_init
    use areaMod   , only : celledge, cellarea                      
    use fileutils , only : getfil
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    type(domain_type),intent(inout) :: domain   ! domain to init
    character(len=*) ,intent(in)    :: filename ! grid filename
!
! !CALLED FROM:
! subroutine surfrd in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!
! !LOCAL VARIABLES:
    integer :: ni,nj               ! size of grid on file
    integer :: ncid,dimid,varid    ! netCDF id's
    integer :: ier                 ! error status
    character(len=256)  :: locfn   ! local file name
    integer :: ret, time_index
    logical :: AREAset             ! true if area read from grid file
    logical :: NSEWset             ! true if lat/lon NSEW read from grid file
    logical :: EDGEset             ! true if EDGE NSEW read from grid file
    integer  :: strt3(3)           ! Start index to read in
    integer  :: cnt3(3)            ! Number of points to read in
    integer  :: strt1, cnt1        ! Start and count to read in for scalar
    real(r8) :: closelat           ! Single-column latitude value
    real(r8) :: closelon           ! Single-column longitude value
    integer  :: closelatidx        ! Single-column latitude index to retrieve
    integer  :: closelonidx        ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_get_grid'     ! subroutine name
!-----------------------------------------------------------------------

    AREAset = .false.
    NSEWset = .false.
    EDGEset = .false.

    if (masterproc) then

       if (filename == ' ') then
          write(6,*) trim(subname),' ERROR: filename must be specified '
          call endrun()
       endif

       call getfil( filename, locfn, 0 )
       call check_ret( nf_open(locfn, 0, ncid), subname )

       if (single_column) then
          ni = lsmlon
          nj = lsmlat
       else
          call check_ret(nf_inq_dimid (ncid, 'lsmlon', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, ni), subname)
          call check_ret(nf_inq_dimid (ncid, 'lsmlat', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, nj), subname)
       endif

    endif

    call mpi_bcast (ni, 1, MPI_INTEGER, 0, mpicom, ier)
    call mpi_bcast (nj, 1, MPI_INTEGER, 0, mpicom, ier)

    call domain_init(domain,ni,nj)

    if (masterproc) then

       domain%edges(:) = spval
       ier = nf_inq_varid (ncid, 'EDGEN', varid)
       if (ier == NF_NOERR) then
          EDGEset = .true.
          call check_ret(nf_inq_varid(ncid, 'EDGEN', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, domain%edges(1)), subname)
          call check_ret(nf_inq_varid(ncid, 'EDGEE', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, domain%edges(2)), subname)
          call check_ret(nf_inq_varid(ncid, 'EDGES', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, domain%edges(3)), subname)
          call check_ret(nf_inq_varid (ncid, 'EDGEW', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, domain%edges(4)), subname)
          if (maxval(domain%edges) > 1.0e35) EDGEset = .false. !read garbage
       endif

       strt3(1)=1
       strt3(2)=1
       strt3(3)=1
       cnt3(1)=domain%ni
       cnt3(2)=domain%nj
       cnt3(3)=1
       strt1=1
       cnt1=domain%nj

       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          strt3(1)=closelonidx
          strt3(2)=closelatidx
          strt1=closelatidx
          cnt1=1
       endif


       ier = nf_inq_varid (ncid, 'LATN', varid)
       if (ier == NF_NOERR) then
          NSEWset = .true.
          call check_ret(nf_inq_varid(ncid, 'LATN', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid,strt3,cnt3, domain%latn), subname)
          call check_ret(nf_inq_varid(ncid, 'LONE', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid,strt3,cnt3, domain%lone), subname)
          call check_ret(nf_inq_varid(ncid, 'LATS', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid,strt3,cnt3, domain%lats), subname)
          call check_ret(nf_inq_varid(ncid, 'LONW', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid,strt3,cnt3, domain%lonw), subname)
       endif

       ier = nf_inq_varid (ncid, 'AREA', varid)
       if (ier == NF_NOERR) then
          AREAset = .true.
          call check_ret(nf_inq_varid(ncid, 'AREA', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid,strt3,cnt3, domain%area), subname)
       endif

       call check_ret(nf_inq_varid(ncid, 'LONGXY' , varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, domain%lonc), subname)
       
       call check_ret(nf_inq_varid(ncid, 'LATIXY', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, domain%latc), subname)
         
! set mask to 1 everywhere by default, override if LANDMASK exists
! if landmask exists, use it to set pftm (for older datasets)
! pftm should be overwritten below for newer datasets
       domain%mask = 1
       ier = nf_inq_varid(ncid, 'LANDMASK', varid)
       if (ier == NF_NOERR) then
          call check_ret(nf_get_vara_int(ncid, varid, strt3, cnt3, domain%mask), subname)
          domain%pftm = domain%mask
          where (domain%mask <= 0)
             domain%pftm = -1
          endwhere
       endif

       ier = nf_inq_varid (ncid, 'PFTDATA_MASK', varid)
       if (ier == NF_NOERR) then
         call check_ret(nf_get_vara_int(ncid, varid, strt3, cnt3, domain%pftm), subname)
       endif
       
!tcx fix, this or a test/abort should be added so overlaps can be computed
!tcx fix, this is demonstrated not bfb in cam bl311 test.
!tcx fix, see also lat_o_local in areaMod.F90
#if (1 == 0)
       ! Check lat limited to -90,90
       if (minval(domain%latc) < -90.0_r8 .or. &
           maxval(domain%latc) >  90.0_r8) then
           write(6,*) trim(subname),' Limiting lat/lon to [-90/90] from ', &
              minval(domain%latc),maxval(domain%latc)
           where (domain%latc < -90.0_r8) domain%latc = -90.0_r8
           where (domain%latc >  90.0_r8) domain%latc =  90.0_r8
       endif
#endif

       ! -------------------------------------------------------------------
       ! Define edges and area of grid cells
       ! -------------------------------------------------------------------

#if (defined SEQ_MCT) || (SEQ_ESMF) || (defined COUP_CSM)
       if (.not.NSEWset) call celledge (domain)
       if (.not.AREAset) call cellarea (domain)
#endif

#if (defined OFFLINE)
       if (.not.NSEWset) then    
          if (.not.EDGEset) then           ! global grid without use of edges
            call celledge (domain)
          else                             ! regional regular grid 
            call celledge (domain, &
                           domain%edges(1), domain%edges(2), &
                           domain%edges(3), domain%edges(4))
          end if	
       endif
       if (.not.AREAset) then
          call cellarea (domain)
       end if	
#endif

       call check_ret(nf_close(ncid), subname)

    end if   ! end of if-masterproc block

    call mpi_bcast (domain%latn , size(domain%latn) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%lats , size(domain%lats) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%lonw , size(domain%lonw) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%lone , size(domain%lone) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%area , size(domain%area) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%latc , size(domain%latc) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%lonc , size(domain%lonc) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (domain%pftm , size(domain%pftm) , MPI_INTEGER, 0, mpicom, ier)
    call mpi_bcast (domain%edges, size(domain%edges), MPI_REAL8  , 0, mpicom, ier)

  end subroutine surfrd_get_grid

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_get_latlon
!
! !INTERFACE:
  subroutine surfrd_get_latlon(latlon,filename,mask,mfilename,pftmflag)
!
! !DESCRIPTION:
! Read the surface dataset grid related information:
!
! !USES:
    use clm_varcon, only : spval
    use domainMod , only : latlon_type, latlon_init
    use areaMod   , only : celledge
    use fileutils , only : getfil
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    type(latlon_type)         ,intent(inout) :: latlon   ! domain to init
    character(len=*)          ,intent(in)    :: filename ! grid filename
    integer,pointer  ,optional               :: mask(:)
    character(len=*) ,optional,intent(in)    :: mfilename ! grid filename
    logical          ,optional,intent(in)    :: pftmflag   ! is mask pft mask?
!
! !CALLED FROM:
! subroutine surfrd in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!
! !LOCAL VARIABLES:
    integer :: ni,nj               ! size of grid on file
    integer :: n                   ! index
    integer :: ncid,dimid,varid    ! netCDF id's
    integer :: ncidm               ! mask file netCDF id's
    integer :: ier                 ! error status
    character(len=256)  :: locfn   ! local file name
    integer :: ret, time_index
    real(r8),pointer :: rdata(:,:) ! temporary data
    logical :: NSEWset             ! true if lat/lon NSEW read from grid file
    logical :: EDGEset             ! true if EDGE NSEW read from grid file
    logical :: lpftmflag           ! is mask a pft mask, local copy
    integer  :: start2(2)          ! Start index to read in
    integer  :: count2(2)          ! Number of points to read in
    real(r8) :: closelat           ! Single-column latitude value
    real(r8) :: closelon           ! Single-column longitude value
    integer  :: closelatidx        ! Single-column latitude index to retrieve
    integer  :: closelonidx        ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_get_latlon'    ! subroutine name
!-----------------------------------------------------------------------

    NSEWset = .false.
    EDGEset = .false.
    lpftmflag = .false.
    ni = 0
    nj = 0

    if (present(pftmflag)) then
       lpftmflag = pftmflag
    endif

    if (masterproc) then

       if (filename == ' ') then
          write(6,*) trim(subname),' ERROR: filename must be specified '
          call endrun()
       endif

       call getfil( filename, locfn, 0 )
       call check_ret( nf_open(locfn, 0, ncid), subname )

       if (single_column) then
          ni = lsmlon
          nj = lsmlat
       else
          ier = nf_inq_dimid (ncid, 'lon', dimid)
          if (ier == NF_NOERR) then
            call check_ret(nf_inq_dimlen(ncid, dimid, ni), subname)
          endif
          ier = nf_inq_dimid (ncid, 'lat', dimid)
          if (ier == NF_NOERR) then
            call check_ret(nf_inq_dimlen(ncid, dimid, nj), subname)
          endif
          ier = nf_inq_dimid (ncid, 'lsmlon', dimid)
          if (ier == NF_NOERR) then
            call check_ret(nf_inq_dimlen(ncid, dimid, ni), subname)
          endif
          ier = nf_inq_dimid (ncid, 'lsmlat', dimid)
          if (ier == NF_NOERR) then
            call check_ret(nf_inq_dimlen(ncid, dimid, nj), subname)
          endif
       endif

       if (ni == 0 .or. nj == 0) then
          write(6,*) trim(subname),' ERROR: ni or nj not set',ni,nj
          call endrun()
       endif

    endif

    call mpi_bcast (ni, 1, MPI_INTEGER, 0, mpicom, ier)
    call mpi_bcast (nj, 1, MPI_INTEGER, 0, mpicom, ier)

    call latlon_init(latlon,ni,nj)
    if (present(mask)) then
       allocate(mask(ni*nj))
       mask = 1
    endif

    if (masterproc) then
       allocate(rdata(ni,nj))

       start2(1)  = 1
       start2(2)  = 1
       count2(1)  = ni
       count2(2)  = nj

       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          start2(1) = closelonidx
          start2(2) = closelatidx
       endif

       call check_ret(nf_inq_varid(ncid, 'LONGXY' , varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
       !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
       latlon%lonc(:) = rdata(:,1)
          
       call check_ret(nf_inq_varid(ncid, 'LATIXY', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
       !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
       latlon%latc(:) = rdata(1,:)

       latlon%edges(:) = spval
       ier = nf_inq_varid (ncid, 'EDGEN', varid)
       if (ier == NF_NOERR) then
          EDGEset = .true.
          call check_ret(nf_inq_varid(ncid, 'EDGEN', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, latlon%edges(1)), subname)
          call check_ret(nf_inq_varid(ncid, 'EDGEE', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, latlon%edges(2)), subname)
          call check_ret(nf_inq_varid(ncid, 'EDGES', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, latlon%edges(3)), subname)
          call check_ret(nf_inq_varid (ncid, 'EDGEW', varid), subname)
          call check_ret(nf_get_var_double(ncid, varid, latlon%edges(4)), subname)
          if (maxval(latlon%edges) > 1.0e35) EDGEset = .false. !read garbage
       endif

       ier = nf_inq_varid (ncid, 'LATN', varid)
       if (ier == NF_NOERR) then
          NSEWset = .true.
          call check_ret(nf_inq_varid(ncid, 'LATN', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
          !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
          latlon%latn(:) = rdata(1,:)

          call check_ret(nf_inq_varid(ncid, 'LONE', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
          !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
          latlon%lone(:) = rdata(:,1)

          call check_ret(nf_inq_varid(ncid, 'LATS', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
          !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
          latlon%lats(:) = rdata(1,:)

          call check_ret(nf_inq_varid(ncid, 'LONW', varid), subname)
          call check_ret(nf_get_vara_double(ncid, varid, start2, count2, rdata), subname)
          !call check_ret(nf_get_var_double(ncid, varid, rdata), subname)
          latlon%lonw(:) = rdata(:,1)
       endif

#if (defined SEQ_MCT) || (defined SEQ_ESMF) || (defined COUP_CSM)
       if (.not.NSEWset) call celledge (latlon)
#endif

#if (defined OFFLINE)
       if (.not.NSEWset) then    
          if (.not.EDGEset) then           ! global grid without use of edges
            call celledge (latlon)
          else                             ! regional regular grid 
            call celledge (latlon, &
                           latlon%edges(1), latlon%edges(2), &
                           latlon%edges(3), latlon%edges(4))
          end if	
       endif
#endif
       if (present(mask)) then
          if (present(mfilename)) then
             if (mfilename == ' ') then
               write(6,*) trim(subname),' ERROR: mfilename must be specified '
               call endrun()
             endif

             call getfil( mfilename, locfn, 0 )
             call check_ret( nf_open(locfn, 0, ncidm), subname )
          else
             ncidm = ncid
          endif

          mask = 1
          ier = nf_inq_varid(ncidm, 'LANDMASK', varid)
          if (ier == NF_NOERR) then
             call check_ret(nf_get_vara_int(ncidm, varid, start2, count2, mask), subname)
          endif

          !--- if this is a pft mask, then modify and look for pftdata_mask array on dataset ---
          if (lpftmflag) then
             do n = 1,ni*nj
                if (mask(n) <= 0) mask(n) = -1
             enddo
             ier = nf_inq_varid (ncidm, 'PFTDATA_MASK', varid)
             if (ier == NF_NOERR) then
                call check_ret(nf_get_vara_int(ncidm, varid, start2, count2, mask), subname)
             endif
          endif

          if (present(mfilename)) then
             call check_ret(nf_close(ncidm), subname)
          endif

       endif

       deallocate(rdata)

!tcx fix, this or a test/abort should be added so overlaps can be computed
!tcx fix, this is demonstrated not bfb in cam bl311 test.
!tcx fix, see also lat_o_local in areaMod.F90
#if (1 == 0)
       ! Check lat limited to -90,90
       if (minval(latlon%latc) < -90.0_r8 .or. &
           maxval(latlon%latc) >  90.0_r8) then
           write(6,*) trim(subname),' Limiting lat/lon to [-90/90] from ', &
              minval(latlon%latc),maxval(latlon%latc)
           where (latlon%latc < -90.0_r8) latlon%latc = -90.0_r8
           where (latlon%latc >  90.0_r8) latlon%latc =  90.0_r8
       endif
#endif

       call check_ret(nf_close(ncid), subname)

    end if   ! end of if-masterproc block

    call mpi_bcast (latlon%latc , size(latlon%latc) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%lonc , size(latlon%lonc) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%lats , size(latlon%lats) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%latn , size(latlon%latn) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%lonw , size(latlon%lonw) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%lone , size(latlon%lone) , MPI_REAL8  , 0, mpicom, ier)
    call mpi_bcast (latlon%edges, size(latlon%edges), MPI_REAL8  , 0, mpicom, ier)
    if (present(mask)) then
       call mpi_bcast(mask       , size(mask)         , MPI_INTEGER, 0, mpicom, ier)
    endif

  end subroutine surfrd_get_latlon

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_get_frac
!
! !INTERFACE:
  subroutine surfrd_get_frac(domain,filename)
!
! !DESCRIPTION:
! Read the landfrac dataset grid related information:
! Assume domain has already been initialized and read
!
! !USES:
    use domainMod , only : domain_type
    use fileutils , only : getfil
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    type(domain_type),intent(inout) :: domain   ! domain to init
    character(len=*) ,intent(in)    :: filename ! grid filename
!
! !CALLED FROM:
! subroutine initialize
!
! !REVISION HISTORY:
! Created by T Craig
!
!EOP
!
! !LOCAL VARIABLES:
    integer :: n                   ! indices
    integer :: ni,nj,ns            ! size of grid on file
    integer :: ncid,dimid,varid    ! netCDF id's
    integer :: ier                 ! error status
    real(r8):: eps = 1.0e-12_r8    ! lat/lon error tolerance
    real(r8),allocatable:: lonc(:),latc(:)  ! local lat/lon
    character(len=256)  :: locfn   ! local file name
    integer :: ret, time_index
    integer  :: strt3(3)           ! Start index to read in
    integer  :: cnt3(3)            ! Number of points to read in
    integer  :: strt1, cnt1        ! Start and count to read in for scalar
    real(r8) :: closelat           ! Single-column latitude value
    real(r8) :: closelon           ! Single-column longitude value
    integer  :: closelatidx        ! Single-column latitude index to retrieve
    integer  :: closelonidx        ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_get_frac'     ! subroutine name
!-----------------------------------------------------------------------

    if (masterproc) then

       if (filename == ' ') then
          write(6,*) trim(subname),' ERROR: filename must be specified '
          call endrun()
       endif

       call getfil( filename, locfn, 0 )
       call check_ret( nf_open(locfn, 0, ncid), subname )

       if (single_column) then
          ni = lsmlon
          nj = lsmlat
       else
          call check_ret(nf_inq_dimid (ncid, 'lsmlon', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, ni), subname)
          call check_ret(nf_inq_dimid (ncid, 'lsmlat', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, nj), subname)
       endif

       ns = ni*nj

       if (domain%ni /= ni .or. domain%nj /= nj .or. domain%ns /= ns) then
          write(6,*) trim(subname),' ERROR: landfrac file mismatch ni,nj',domain%ni,ni,domain%nj,nj,domain%ns,ns
          call endrun()
       endif

       allocate(latc(ni*nj),lonc(ni*nj))

       strt3(1)=1
       strt3(2)=1
       strt3(3)=1
       cnt3(1)=domain%ni
       cnt3(2)=domain%nj
       cnt3(3)=1
       strt1=1
       cnt1=domain%nj

       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          strt3(1)=closelonidx
          strt3(2)=closelatidx
          strt1=closelatidx
          cnt1=1
       endif

       call check_ret(nf_inq_varid(ncid, 'LONGXY' , varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, lonc), subname)
          
       call check_ret(nf_inq_varid(ncid, 'LATIXY', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, latc), subname)

       do n = 1,ns
          if (abs(latc(n)-domain%latc(n)) > eps .or. &
               abs(lonc(n)-domain%lonc(n)) > eps) then
             write(6,*) trim(subname),' ERROR: landfrac file mismatch lat,lon',latc(n),domain%latc(n),lonc(n),domain%lonc(n),eps
             call endrun()
          endif
       enddo
       
       call check_ret(nf_inq_varid(ncid, 'LANDMASK', varid), subname)
       call check_ret(nf_get_vara_int(ncid, varid, strt3, cnt3, domain%mask), subname)
       
       call check_ret(nf_inq_varid(ncid, 'LANDFRAC', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, domain%frac), subname)
       

       deallocate(latc,lonc)

       call check_ret(nf_close(ncid), subname)

    end if   ! end of if-masterproc block

    call mpi_bcast (domain%mask , size(domain%mask) , MPI_INTEGER, 0, mpicom, ier)
    call mpi_bcast (domain%frac , size(domain%frac) , MPI_REAL8  , 0, mpicom, ier)

  end subroutine surfrd_get_frac

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_get_topo
!
! !INTERFACE:
  subroutine surfrd_get_topo(domain,filename)
!
! !DESCRIPTION:
! Read the topo dataset grid related information:
! Assume domain has already been initialized and read
!
! !USES:
    use domainMod , only : domain_type
    use fileutils , only : getfil
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    type(domain_type),intent(inout) :: domain   ! domain to init
    character(len=*) ,intent(in)    :: filename ! grid filename
!
! !CALLED FROM:
! subroutine initialize
!
! !REVISION HISTORY:
! Created by T Craig
!
!EOP
!
! !LOCAL VARIABLES:
    integer :: n                   ! indices
    integer :: ni,nj,ns            ! size of grid on file
    integer :: ncid,dimid,varid    ! netCDF id's
    integer :: ier                 ! error status
    real(r8):: eps = 1.0e-12_r8    ! lat/lon error tolerance
    real(r8),allocatable:: lonc(:),latc(:)  ! local lat/lon
    character(len=256)  :: locfn   ! local file name
    integer :: ret, time_index
    integer  :: strt3(3)           ! Start index to read in
    integer  :: cnt3(3)            ! Number of points to read in
    integer  :: strt1, cnt1        ! Start and count to read in for scalar
    real(r8) :: closelat           ! Single-column latitude value
    real(r8) :: closelon           ! Single-column longitude value
    integer  :: closelatidx        ! Single-column latitude index to retrieve
    integer  :: closelonidx        ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_get_topo'     ! subroutine name
!-----------------------------------------------------------------------

    if (masterproc) then

       if (filename == ' ') then
          write(6,*) trim(subname),' ERROR: filename must be specified '
          call endrun()
       endif

       call getfil( filename, locfn, 0 )
       call check_ret( nf_open(locfn, 0, ncid), subname )

       if (single_column) then
          ni = lsmlon
          nj = lsmlat
       else
          call check_ret(nf_inq_dimid (ncid, 'lsmlon', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, ni), subname)
          call check_ret(nf_inq_dimid (ncid, 'lsmlat', dimid), subname)
          call check_ret(nf_inq_dimlen(ncid, dimid, nj), subname)
       endif

       ns = ni*nj

       if (domain%ni /= ni .or. domain%nj /= nj .or. domain%ns /= ns) then
          write(6,*) trim(subname),' ERROR: topo file mismatch ni,nj',domain%ni,ni,domain%nj,nj,domain%ns,ns
          call endrun()
       endif

       allocate(latc(ns),lonc(ns))

       strt3(1)=1
       strt3(2)=1
       strt3(3)=1
       cnt3(1)=domain%ni
       cnt3(2)=domain%nj
       cnt3(3)=1
       strt1=1
       cnt1=domain%nj

       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          strt3(1)=closelonidx
          strt3(2)=closelatidx
          strt1=closelatidx
          cnt1=1
       endif

       call check_ret(nf_inq_varid(ncid, 'LONGXY' , varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, lonc), subname)

       call check_ret(nf_inq_varid(ncid, 'LATIXY', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, latc), subname)

       do n = 1,ns
          if (abs(latc(n)-domain%latc(n)) > eps .or. &
              abs(lonc(n)-domain%lonc(n)) > eps) then
             write(6,*) trim(subname),' ERROR: topo file mismatch lat,lon',latc(n),domain%latc(n),lonc(n),domain%lonc(n),eps
             call endrun()
          endif
       enddo

       call check_ret(nf_inq_varid(ncid, 'TOPO', varid), subname)
       call check_ret(nf_get_vara_double(ncid, varid, strt3, cnt3, domain%topo), subname)

       deallocate(latc,lonc)

       call check_ret(nf_close(ncid), subname)

    end if   ! end of if-masterproc block

    call mpi_bcast (domain%topo , size(domain%topo) , MPI_REAL8  , 0, mpicom, ier)

  end subroutine surfrd_get_topo

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_wtxy_special 
!
! !INTERFACE:
!  subroutine surfrd_wtxy_special(ncid, pctspec, vegxy, wtxy, domain)
  subroutine surfrd_wtxy_special(ncid, domain)
!
! !DESCRIPTION:
! Determine weight with respect to gridcell of all special "pfts" as well
! as soil color and percent sand and clay
!
! !USES:
    use pftvarcon   , only : noveg
    use domainMod   , only : domain_type
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    integer , intent(in)    :: ncid      ! netcdf file id 
!    real(r8), intent(inout) :: pctspec(:)! percent wrt gcell special lunits
!    integer , intent(inout) :: vegxy(:,:)  ! PFT
!    real(r8), intent(inout) :: wtxy(:,:) ! subgrid weights
    type(domain_type),intent(in) :: domain ! domain associated with wtxy
!
! !CALLED FROM:
! subroutine surfrd in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein, Sam Levis and Gordon Bonan
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: n,nl,ns                    ! indices
    integer  :: begg,endg                  ! gcell beg/end
    integer  :: dimid,varid                ! netCDF id's
    integer  :: ret, time_index
    real(r8) :: nlevsoidata(nlevsoi)
    logical  :: found                      ! temporary for error check
    integer  :: nindx                      ! temporary for error check
    integer  :: ier                        ! error status
    real(r8),pointer :: pctgla(:)      ! percent of grid cell is glacier
    real(r8),pointer :: pctlak(:)      ! percent of grid cell is lake
    real(r8),pointer :: pctwet(:)      ! percent of grid cell is wetland
    real(r8),pointer :: pcturb(:)      ! percent of grid cell is urbanized
    integer  :: strt3(3)               ! Start index to read in
    integer  :: cnt3(3)                ! Number of points to read in
    real(r8) :: closelat               ! Single-column latitude value
    real(r8) :: closelon               ! Single-column longitude value
    integer  :: closelatidx            ! Single-column latitude index to retrieve
    integer  :: closelonidx            ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_wtxy_special'  ! subroutine name
!!-----------------------------------------------------------------------

    ns = domain%ns
    call get_proc_bounds(begg,endg)

    allocate(pctgla(begg:endg),pctlak(begg:endg))
    allocate(pctwet(begg:endg),pcturb(begg:endg))

    if (masterproc) then
       call check_dim(ncid, 'nlevsoi', nlevsoi)
    end if   ! end of if-masterproc

    ! Obtain non-grid surface properties of surface dataset other than percent pft

    strt3(1)=1
    strt3(2)=1
    strt3(3)=1
    cnt3(1)=domain%ni
    cnt3(2)=domain%nj
    cnt3(3)=1

    if (single_column) then
       call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
       strt3(1)=closelonidx
       strt3(2)=closelatidx
    endif

    call ncd_iolocal(ncid, 'PCT_WETLAND', 'read', pctwet, begg, endg, gsMap_lnd_gdc2glo, perm_lnd_gdc2glo, strt3, cnt3 )
    call ncd_iolocal(ncid, 'PCT_LAKE'   , 'read', pctlak, begg, endg, gsMap_lnd_gdc2glo, perm_lnd_gdc2glo, strt3, cnt3 )
    call ncd_iolocal(ncid, 'PCT_GLACIER', 'read', pctgla, begg, endg, gsMap_lnd_gdc2glo, perm_lnd_gdc2glo, strt3, cnt3 )
    call ncd_iolocal(ncid, 'PCT_URBAN'  , 'read', pcturb, begg, endg, gsMap_lnd_gdc2glo, perm_lnd_gdc2glo, strt3, cnt3 )
    pctspec = pctwet + pctlak + pctgla + pcturb

    ! Error check: glacier, lake, wetland, urban sum must be less than 100

    found = .false.
    do nl = begg,endg
       if (pctspec(nl) > 100._r8+1.e-04_r8) then
          found = .true.
          nindx = nl
          exit
       end if
       if (found) exit
    end do
    if ( found ) then
       write(6,*)'surfrd error: PFT cover>100 for nl=',nindx
       call endrun()
    end if

    ! Error check that urban parameterization is not yet finished

    found = .false.
    do nl = begg,endg
       if (pcturb(nl) /= 0._r8) then
          found = .true.
          nindx = nl
          exit
       end if
       if (found) exit
    end do
    if ( found ) then
       write (6,*)'surfrd error: urban parameterization not implemented at nl= ',nindx,pcturb(nl)
       call endrun()
    end if

    ! Determine veg and wtxy for special landunits

    do nl = begg,endg
       vegxy(nl,npatch_urban)  = noveg
       wtxy(nl,npatch_urban)   = pcturb(nl)/100._r8

       vegxy(nl,npatch_lake)   = noveg
       wtxy(nl,npatch_lake)    = pctlak(nl)/100._r8

       vegxy(nl,npatch_wet)    = noveg
       wtxy(nl,npatch_wet)     = pctwet(nl)/100._r8

       vegxy(nl,npatch_glacier)= noveg
       wtxy(nl,npatch_glacier) = pctgla(nl)/100._r8
    end do

   deallocate(pctgla,pctlak,pctwet,pcturb)

  end subroutine surfrd_wtxy_special

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_wtxy_veg_rank
!
! !INTERFACE:
!  subroutine surfrd_wtxy_veg_rank(ncid, pctspec, vegxy, wtxy, domain)
  subroutine surfrd_wtxy_veg_rank(ncid, domain)
!
! !DESCRIPTION:
! Determine wtxy and veg arrays for non-dynamic landuse mode
!
! !USES:
    use clm_varctl, only : create_crop_landunit
    use pftvarcon   , only : crop, noveg
    use domainMod   , only : domain_type
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    integer , intent(in)    :: ncid        ! netcdf file id 
!    real(r8), intent(in)    :: pctspec(:)  ! percent wrt gcell of spec lunits
!    integer , intent(inout) :: vegxy(:,:)    ! PFT
!    real(r8), intent(inout) :: wtxy(:,:)   ! subgrid weights
    type(domain_type),intent(in) :: domain ! domain associated with wtxy
!
! !CALLED FROM:
! subroutine surfrd in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein, Sam Levis and Gordon Bonan
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: k,m,k1,k2,n,nl,ns               ! indices
    integer  :: begg,endg                       ! beg/end gcell index
    integer  :: dimid,varid                     ! netCDF id's
    integer  :: start(3),count(3)               ! netCDF start/count arrays
    integer  :: cropcount                       ! temporary counter
    real(r8),allocatable :: sumvec(:)           ! temporary vector sum
    logical  :: found                           ! temporary for error check
    integer  :: nindx                           ! temporary for error check
    integer  :: miss = 99999                    ! missing data indicator
    real(r8) :: wst(0:numpft)                   ! as pft at specific i, j
    integer ,allocatable :: wsti(:)             ! ranked indices largest wst values
    real(r8) :: wst_sum                         ! sum of %pft
    real(r8) :: sumpct                          ! sum of %pft over maxpatch_pft
    real(r8) :: diff                            ! the difference (wst_sum - sumpct)
    real(r8) :: rmax                            ! maximum patch cover
    integer ,allocatable :: pft(:,:)            ! PFT
    integer ,allocatable :: cft(:,:)            ! CFT
    real(r8),allocatable :: pctcft_lunit(:,:)   ! % of crop lunit area for CFTs
    real(r8),allocatable :: pctpft_lunit(:,:)   ! % of vegetated lunit area PFTs
    integer  :: ier                             ! error status
    real(r8),allocatable :: pctpft(:,:)         ! percent of vegetated gridcell area for PFTs
    real(r8),pointer :: arrayl(:)               ! local array
    integer ,pointer :: irrayg(:)               ! global array
    integer  :: ret, time_index
    real(r8),allocatable :: rmaxpatchdata(:)
    integer ,allocatable :: imaxpatchdata(:)
    real(r8) :: numpftp1data(0:numpft)         
    integer  :: strt3(3)                        ! Start index to read in
    integer  :: cnt3(3)                         ! Number of points to read in
    real(r8) :: closelat                        ! Single-column latitude value
    real(r8) :: closelon                        ! Single-column longitude value
    integer  :: closelatidx                     ! Single-column latitude index to retrieve
    integer  :: closelonidx                     ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_wtxy_veg_rank'  ! subroutine name
!-----------------------------------------------------------------------

    ns = domain%ns
    call get_proc_bounds(begg,endg)

    allocate(sumvec(begg:endg))
    allocate(cft(begg:endg,maxpatch_cft))
    allocate(pft(begg:endg,maxpatch_pft))
    allocate(pctcft_lunit(begg:endg,maxpatch_cft))
    allocate(pctpft_lunit(begg:endg,maxpatch_pft))
    allocate(pctpft(begg:endg,0:numpft))
    allocate(wsti(maxpatch_pft))

    if (masterproc) then
       call check_dim(ncid, 'lsmpft', numpft+1)
    end if
    allocate(arrayl(begg:endg))
    do n = 0,numpft
       start(1) = 1
       count(1) = domain%ni
       start(2) = 1
       count(2) = domain%nj
       start(3) = n+1	 ! dataset is 1:numpft+1, not 0:numpft
       count(3) = 1
       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          start(1)=closelonidx
          start(2)=closelatidx
       endif
       call ncd_iolocal(ncid, 'PCT_PFT', 'read', arrayl, begg, endg, gsMap_lnd_gdc2glo, perm_lnd_gdc2glo, start, count)
       pctpft(begg:endg,n) = arrayl(begg:endg)
    enddo
    deallocate(arrayl)

    ! 1. pctpft must go back to %vegetated landunit instead of %gridcell
    ! 2. pctpft bare = 100 when landmask = 1 and 100% special landunit
    ! NB: (1) and (2) do not apply to crops.
    ! For now keep all cfts (< 4 anyway) instead of 4 most dominant cfts

    do nl = begg,endg

       cft(nl,:) = 0
       pctcft_lunit(nl,:) = 0._r8
       cropcount = 0

       if (pctspec(nl) < 100._r8) then

          do m = 0, numpft
             if (create_crop_landunit) then
                ! Separate crop landunit is to be created

                if (crop(m) == 1._r8 .and. pctpft(nl,m) > 0._r8) then
                   cropcount = cropcount + 1
                   if (cropcount > maxpatch_cft) then
                      write(6,*) 'ERROR surfrdMod: cropcount>maxpatch_cft'
                      call endrun()
                   end if
                   cft(nl,cropcount) = m
                   pctcft_lunit(nl,cropcount) = pctpft(nl,m) * 100._r8/(100._r8-pctspec(nl))
                   pctpft(nl,m) = 0.0_r8
                else if (crop(m) == 0._r8) then
                   pctpft(nl,m) = pctpft(nl,m) * 100._r8/(100._r8-pctspec(nl))
                end if

             else
                ! Separate crop landunit is not created

                pctpft(nl,m) = pctpft(nl,m) * 100._r8/(100._r8-pctspec(nl))
             end if
          end do

       else if (pctspec(nl) == 100._r8) then

          pctpft(nl,0)        = 100._r8
          pctpft(nl,1:numpft) =   0._r8

       else

          write(6,*)subname, 'error: pcturb+pctgla+pctlak+pctwet = ',pctspec(nl), &
               ' must be less than or equal to 100'
          call endrun()

       end if
    end do

    ! Find pft and pct arrays 
    ! Save percent cover by PFT [wst] and total percent cover [wst_sum]

    do nl = begg,endg

       wst_sum = 0._r8
       sumpct = 0
       do m = 0, numpft
          wst(m) = pctpft(nl,m)
          wst_sum = wst_sum + pctpft(nl,m)
       end do

       if (domain%pftm(nl) >= 0) then

          ! Rank [wst] in ascendg order to obtain the top [maxpatch_pft] PFTs

          call surfrd_mkrank (numpft, wst, miss, wsti, maxpatch_pft)

          ! Fill in [pft] and [pctpft] with data for top [maxpatch_pft] PFTs.
          ! If land model grid cell is ocean, set to no PFTs.
          ! If land model grid cell is land then:
          !  1. If [pctlnd_o] = 0, there is no PFT data from the input grid.
          !     Since need land data, use bare ground.
          !  2. If [pctlnd_o] > 0, there is PFT data from the input grid but:
          !     a. use the chosen PFT so long as it is not a missing value
          !     b. missing value means no more PFTs with cover > 0
                   
          do m = 1, maxpatch_pft
             if (wsti(m) /=  miss) then
                pft(nl,m) = wsti(m)
                pctpft_lunit(nl,m) = wst(wsti(m))
             else
                pft(nl,m) = noveg
                pctpft_lunit(nl,m) = 0._r8
             end if
             sumpct = sumpct + pctpft_lunit(nl,m)
          end do

       else                               ! model grid wants ocean

          do m = 1, maxpatch_pft
             pft(nl,m) = 0
             pctpft_lunit(nl,m) = 0._r8
          end do

       end if

       ! Correct for the case of more than [maxpatch_pft] PFTs present
                
       if (sumpct < wst_sum) then
          diff  = wst_sum - sumpct
          sumpct = 0._r8
          do m = 1, maxpatch_pft
             pctpft_lunit(nl,m) = pctpft_lunit(nl,m) + diff/maxpatch_pft
             sumpct = sumpct + pctpft_lunit(nl,m)
          end do
       end if

       ! Error check: make sure have a valid PFT

       do m = 1,maxpatch_pft
          if (pft(nl,m) < 0 .or. pft(nl,m) > numpft) then
             write (6,*)'surfrd error: invalid PFT at gridcell nl=',nl,pft(nl,m)
             call endrun()
          end if
       end do

       ! As done in mksrfdatMod.F90 for other percentages, truncate pctpft to
       ! ensure that weight relative to landunit is not nonzero
       ! (i.e. a very small number such as 1e-16) where it really should be zero
       ! The following if-block is here to preserve roundoff level differences
       ! between the call to surfrd_wtxy_veg_all and surfrd_wtxy_veg_rank

       if (maxpatch_pft < numpft+1) then
          do m=1,maxpatch_pft
             pctpft_lunit(nl,m) = float(nint(pctpft_lunit(nl,m)))
          end do
          do m=1,maxpatch_cft
             pctcft_lunit(nl,m) = float(nint(pctcft_lunit(nl,m)))
          end do
       end if
                   
       ! Make sure sum of PFT cover == 100 for land points. If not,
       ! subtract excess from most dominant PFT.

       rmax = -9999._r8
       k1 = -9999
       k2 = -9999
       sumpct = 0._r8
       do m = 1, maxpatch_pft
          sumpct = sumpct + pctpft_lunit(nl,m)
          if (pctpft_lunit(nl,m) > rmax) then
             k1 = m
             rmax = pctpft_lunit(nl,m)
          end if
       end do
       do m = 1, maxpatch_cft
          sumpct = sumpct + pctcft_lunit(nl,m)
          if (pctcft_lunit(nl,m) > rmax) then
             k2 = m
             rmax = pctcft_lunit(nl,m)
          end if
       end do
       if (k1 == -9999 .and. k2 == -9999) then
          write(6,*)'surfrd error: largest PFT patch not found'
          call endrun()
       else if (domain%pftm(nl) >= 0) then
          if (sumpct < 95 .or. sumpct > 105._r8) then
             write(6,*)'surfrd error: sum of PFT cover =',sumpct,' at nl=',nl
             call endrun()
          else if (sumpct /= 100._r8 .and. k2 /= -9999) then
             pctcft_lunit(nl,k2) = pctcft_lunit(nl,k2) - (sumpct-100._r8)
          else if (sumpct /= 100._r8) then
             pctpft_lunit(nl,k1) = pctpft_lunit(nl,k1) - (sumpct-100._r8)
          end if
       end if

       ! Error check: make sure PFTs sum to 100% cover

       sumpct = 0._r8
       do m = 1, maxpatch_pft
          sumpct = sumpct + pctpft_lunit(nl,m)
       end do
       do m = 1, maxpatch_cft
          sumpct = sumpct + pctcft_lunit(nl,m)
       end do
       if (domain%pftm(nl) >= 0) then
          if (abs(sumpct - 100._r8) > 0.000001_r8) then
             write(6,*)'surfrdMod error: sum(pct) over maxpatch_pft is not = 100.'
             write(6,*)sumpct, nl
             call endrun()
          end if
          if (sumpct < -0.000001_r8) then
             write(6,*)'surfrdMod error: sum(pct) over maxpatch_pft is < 0.'
             write(6,*)sumpct, nl
             call endrun()
          end if
       end if

    end do   ! end of latitude loop

    ! Determine array [veg], which sets the PFT type for each of the [maxpatch]
    ! patches and array [wtxy], which sets the relative abundance of the PFT.
    ! Fill in PFTs for vegetated portion of grid cell. Fractional areas for
    ! these points [pctpft] pertain to "vegetated" area not to total grid area.
    ! So need to adjust them for fraction of grid that is vegetated.
    ! Next, fill in urban, lake, wetland, and glacier patches.

    do nl = begg,endg
       if (domain%pftm(nl) >= 0) then

          ! Naturally vegetated landunit

          do m = 1, maxpatch_pft
             vegxy(nl,m)  = pft(nl,m)
             wtxy(nl,m) = pctpft_lunit(nl,m) * (100._r8-pctspec(nl))/10000._r8
#if (defined CN)
             ! the following test prevents the assignment of temperate deciduous
             ! vegetation types in the tropics
             ! 1. broadleaf deciduous temperate tree -> broadleaf deciduous tropical tree

             if (vegxy(nl,m) == 7 .and. abs(domain%latc(nl)) < 23.5_r8) vegxy(nl,m) = 6

             ! 2. broadleaf deciduous temperate shrub -> broadleaf deciduous tropical tree
             ! this reassignment from shrub to tree is necessary because there is currently no
             ! tropical deciduous broadleaf shrub type defined.

              if (vegxy(nl,m) == 10 .and. abs(domain%latc(nl)) < 23.5_r8) vegxy(nl,m) = 6
#endif
          end do

          ! Crop landunit

          if (create_crop_landunit) then
             do m = 1,maxpatch_cft
                vegxy(nl,npatch_glacier+m)  = cft(nl,m)
                wtxy(nl,npatch_glacier+m) = pctcft_lunit(nl,m) * (100._r8-pctspec(nl))/10000._r8
             end do
          end if

       end if
    end do

    ! Error check

    found = .false.
    sumvec(:) = abs(sum(wtxy,dim=2)-1._r8)
    do nl = begg,endg
       if (sumvec(nl) > 1.e-06_r8 .and. domain%pftm(nl)>=0) then
          found = .true.
          nindx = nl
          exit
       endif
    end do
    if ( found ) then
       write (6,*)'surfrd error: WTXY > 1 occurs at nl= ',nindx; call endrun()
    end if

    deallocate(sumvec,cft,pft)
    deallocate(pctcft_lunit,pctpft_lunit,pctpft)
    deallocate(wsti)

  end subroutine surfrd_wtxy_veg_rank

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_wtxy_veg_all
!
! !INTERFACE:
!  subroutine surfrd_wtxy_veg_all(ncid, pctspec, vegxy, wtxy, domain)
  subroutine surfrd_wtxy_veg_all(ncid, domain)
!
! !DESCRIPTION:
! Determine wtxy and veg arrays for non-dynamic landuse mode
!
! !USES:
    use domainMod   , only : domain_type
!
! !ARGUMENTS:
    implicit none
    include 'netcdf.inc'
    integer , intent(in)    :: ncid       ! netcdf file id 
!    real(r8), intent(in)    :: pctspec(:) ! percent wrt gcell of spec lunits
!    integer , intent(inout) :: vegxy(:,:)   ! PFT
!    real(r8), intent(inout) :: wtxy(:,:)  ! subgrid weights
    type(domain_type),intent(in) :: domain ! domain associated with wtxy
!
! !CALLED FROM:
! subroutine surfrd in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein, Sam Levis and Gordon Bonan
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: m,mp7,mp8,mp11,n,nl,ns         ! indices
    integer  :: begg,endg                      ! beg/end gcell index
    integer  :: dimid,varid                    ! netCDF id's
    integer  :: start(3),count(3)              ! netcdf start/count arrays
    integer  :: ier                            ! error status	
    real(r8) :: sumpct                         ! sum of %pft over maxpatch_pft
    real(r8),allocatable :: pctpft(:,:)        ! percent of vegetated gridcell area for PFTs
    real(r8),pointer :: arrayl(:)              ! local array
    integer  :: ret, time_index
    real(r8) :: numpftp1data(0:numpft)         
    integer  :: strt3(3)                        ! Start index to read in
    integer  :: cnt3(3)                         ! Number of points to read in
    real(r8) :: closelat                        ! Single-column latitude value
    real(r8) :: closelon                        ! Single-column longitude value
    integer  :: closelatidx                     ! Single-column latitude index to retrieve
    integer  :: closelonidx                     ! Single-column longitude index to retrieve
    character(len=32) :: subname = 'surfrd_wtxy_veg_all'  ! subroutine name
!-----------------------------------------------------------------------

    ns = domain%ns
    call get_proc_bounds(begg,endg)
    allocate(pctpft(begg:endg,0:numpft))

    if (masterproc) then
       call check_dim(ncid, 'lsmpft', numpft+1)
    endif

    allocate(arrayl(begg:endg))
    do n = 0,numpft
       start(1) = 1
       count(1) = domain%ni
       start(2) = 1
       count(2) = domain%nj
       start(3) = n+1         ! dataset is 1:numpft+1, not 0:numpft
       count(3) = 1
       if (single_column) then
          call scam_setlatlonidx(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx)
          start(1)=closelonidx
          start(2)=closelatidx
       endif
       call ncd_iolocal(ncid, 'PCT_PFT', 'read', arrayl, begg, endg, gsMap_lnd_gdc2glo, &
                        perm_lnd_gdc2glo, start, count)
       pctpft(begg:endg,n) = arrayl(begg:endg)
    enddo
    deallocate(arrayl)

    do nl = begg,endg
       if (domain%pftm(nl) >= 0) then

          ! Error check: make sure PFTs sum to 100% cover for vegetated landunit 
          ! (convert pctpft from percent with respect to gridcel to percent with 
          ! respect to vegetated landunit)

          if (pctspec(nl) < 100._r8) then
             sumpct = 0._r8
             do m = 0,numpft
                sumpct = sumpct + pctpft(nl,m) * 100._r8/(100._r8-pctspec(nl))
             end do
             if (abs(sumpct - 100._r8) > 0.1e-4_r8) then
                write(6,*)'surfrdMod error: sum(pct) over numpft+1 is not = 100.'
                write(6,*) sumpct, sumpct-100._r8, nl
                call endrun()
             end if
             if (sumpct < -0.000001_r8) then
                write(6,*)'surfrdMod error: sum(pct) over numpft+1 is < 0.'
                write(6,*) sumpct, nl
                call endrun()
             end if
          end if

          ! Set weight of each pft wrt gridcell (note that maxpatch_pft = numpft+1 here)

          do m = 1,numpft+1
             vegxy(nl,m)  = m-1
             wtxy(nl,m) = pctpft(nl,m-1) / 100._r8
          end do

#if (defined CN)
          ! the following test prevents the assignment of temperate deciduous
          ! vegetation types in the tropics
          ! 1. broadleaf deciduous temperate tree (type7) -> broadleaf deciduous tropical tree (type6)
          ! N.B. the veg and wtxy arrays start at 1, so index 1 corresponds to
          ! veg type 0.  So in this case I want to trap on veg types 7 and 10, 
          ! which are indices 8 and 11. Moving to vegtype6, or index 7
          mp7 = 7
          mp8 = 8
          mp11 = 11
          if (abs(domain%latc(nl)) < 23.5_r8 .and. wtxy(nl,mp8) > 0._r8) then
             if (masterproc) then
                write(6,*)'surfrdMod warning: reassigning temperate tree -> tropical tree'
                write(6,*)'nl,lat,veg7wt,veg6wt,type'
                write(6,*) nl,domain%latc(nl),wtxy(nl,mp8),wtxy(nl,mp7),vegxy(nl,mp8)
             end if
             wtxy(nl,mp7) = wtxy(nl,mp7) + wtxy(nl,mp8)
             wtxy(nl,mp8) = 0._r8
             if (masterproc) then
                write(6,*) nl,domain%latc(nl),wtxy(nl,mp8),wtxy(nl,mp7)
             end if
          end if

          ! 2. broadleaf deciduous temperate shrub (type10) -> broadleaf deciduous tropical tree (type6)
          ! this reassignment from shrub to tree is necessary because there is currently no
          ! tropical deciduous broadleaf shrub type defined.

          if (abs(domain%latc(nl)) < 23.5_r8 .and. wtxy(nl,mp11) > 0._r8) then
             if (masterproc) then
                write(6,*)'surfrdMod warning: reassigning temperate shrub -> tropical tree'
                write(6,*)'nl,lat,veg10wt,veg6wt,type'
                write(6,*) nl,domain%latc(nl),wtxy(nl,mp11),wtxy(nl,mp7),vegxy(nl,mp11)
             end if
             wtxy(nl,mp7) = wtxy(nl,mp7) + wtxy(nl,mp11)
             wtxy(nl,mp11) = 0._r8
             if (masterproc) then
                write(6,*) nl,domain%latc(nl),wtxy(nl,mp11),wtxy(nl,mp7)
             end if
          end if
#endif
       end if
    end do

    deallocate(pctpft)

  end subroutine surfrd_wtxy_veg_all

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: surfrd_wtxy_veg_dgvm
!
! !INTERFACE:
!  subroutine surfrd_wtxy_veg_dgvm(pctspec, vegxy, wtxy, domain)
  subroutine surfrd_wtxy_veg_dgvm(domain)
!
! !DESCRIPTION:
! Determine wtxy and vegxy for DGVM mode.
!
! !USES:
    use pftvarcon   , only : crop, noveg
    use domainMod   , only : domain_type
!
! !ARGUMENTS:
    implicit none
!    real(r8), intent(in)    :: pctspec(:) ! percent gridcell of special landunits
!    integer , intent(inout) :: vegxy(:,:)   ! PFT
!    real(r8), intent(inout) :: wtxy(:,:)  ! subgrid weights
    type(domain_type),intent(in) :: domain ! domain associated with wtxy
!
! !CALLED FROM:
! subroutine surfrd in this module
!
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/04
!
!EOP
!
! !LOCAL VARIABLES:
    integer :: m,nl         ! indices
    integer  :: begg,endg   ! beg/end gcell index
!-----------------------------------------------------------------------

    call get_proc_bounds(begg,endg)
    do nl = begg,endg
       do m = 1, maxpatch_pft
          vegxy(nl,m)  = noveg 
          wtxy(nl,m) = 1.0_r8/maxpatch_pft * (100._r8-pctspec(nl))/100._r8
       end do
    end do

  end subroutine surfrd_wtxy_veg_dgvm
   
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: surfrd_mkrank
!
! !INTERFACE:
  subroutine surfrd_mkrank (n, a, miss, iv, num)
!
! !DESCRIPTION:
! Return indices of largest [num] values in array [a]
!
! !USES:
    use shr_kind_mod, only: r8 => shr_kind_r8
    use abortutils, only : endrun
!
! !ARGUMENTS:
    implicit none
    integer , intent(in) :: n        ! array length
    real(r8), intent(in) :: a(0:n)   ! array to be ranked
    integer , intent(in) :: miss     ! missing data value
    integer , intent(in) :: num      ! number of largest values requested
    integer , intent(out):: iv(num)  ! index to [num] largest values in array [a]
!
! !CALLED FROM:
! ! subroutine surfrd_wtxy_veg_rank in this module
!
! !REVISION HISTORY:
! Author: Gordon Bonan
!
!EOP
!
! !LOCAL VARIABLES:
    real(r8) :: a_max       ! maximum value in array
    integer  :: i           ! array index
    real(r8) :: delmax      ! tolerance for finding if larger value
    integer  :: m           ! do loop index
    integer  :: k           ! do loop index
    logical  :: exclude     ! true if data value has already been chosen
!-----------------------------------------------------------------------

    delmax = 1.e-06_r8

    ! Find index of largest non-zero number
    
    iv(1) = miss
    a_max = -9999._r8

    do i = 0, n
       if (a(i)>0._r8 .and. (a(i)-a_max)>delmax) then
          a_max = a(i)
          iv(1)  = i
       end if
    end do

    ! iv(1) = miss indicates no values > 0. this is an error

    if (iv(1) == miss) then
       write (6,*) 'surfrd_mkrank error: iv(1) = missing'
       call endrun
    end if

    ! Find indices of the next [num]-1 largest non-zero number.
    ! iv(m) = miss if there are no more values > 0

    do m = 2, num
       iv(m) = miss
       a_max = -9999._r8
       do i = 0, n

          ! exclude if data value has already been chosen

          exclude = .false.
          do k = 1, m-1
             if (i == iv(k)) exclude = .true.
          end do

          ! if not already chosen, see if it is the largest of
          ! the remaining values

          if (.not. exclude) then
             if (a(i)>0._r8 .and. (a(i)-a_max)>delmax) then
                a_max = a(i)
                iv(m)  = i
             end if
          end if
       end do
    end do

  end subroutine surfrd_mkrank

end module surfrdMod
